
cnvphmmcopy <- function (correctOutput, segmentOutput, chr = c(1:22, "X", "Y"), sample_name= 'cell(s)', ...) 
{
    if (is.null(segmentOutput$segs)) {
        warning("Processed segments now found, automatically processing")
        segmentOutput$segs <- processSegments(segments$segs, 
            space(correctOutput), start(correctOutput), end(correctOutput), 
            correctOutput$copy)
    }
	
	stat_offset = 0
	ylim1=c(-1, 7)
	h_line=c( (ylim1[1]+1):(ylim1[2]-1) )
	chr=as.character(chr)
    segs <- segmentOutput$segs
    correctOutput$state <- segmentOutput$state
    cols <- stateCols()
    # range <- quantile(correctOutput$copy, na.rm = TRUE, prob = c(0.01, 0.99))

	
    a <- correctOutput[chr]
	# a <- correctOutput[which(correctOutput$space %in% chr),]
    # b <- segs[which(segs$chr == chr[1]), ]
	b <- segs[which(segs$chr %in% chr), ]
	data_path=Sys.getenv('data_path')
	chr_len=paste(data_path, 'b37_chr_ln_y.RData', sep='/')
	load(chr_len)
	target_chr_len=b37_chr_ln[which(b37_chr_ln$chr %in% chr),]
	xlim=c(1, sum(as.numeric(target_chr_len$len)))
	ranges=as.data.frame(a$ranges)
	a$x=ranges$start

	pdf("1.cnv.hmmcopy.pdf", width=18, height=6)
	par(pch=18, lwd = .2, ann=F, xaxs='i', yaxs='i', fig=c(0, 1, 0.4, 1))
	
	plot(NA, type='n', xlim=xlim, ylim=ylim1, cex.axis=.5)
	abline(h=h_line, lty=3, col='gray')
	
	
	# plot(a$x[which(a$space==chr[1])], a$copy[which(a$space==chr[1])], col = cols[as.numeric(as.character(a$state))], xlim=xlim, ylim = range, xlab="Position on the genome(bp)", ylab="Copy number states", ...)
	main_title=paste("CNVs of sample ", sample_name)
	title(main=main_title, xlab="Position on the genome(bp)", ylab="Copy number")
	
		
	for (i in 1:nrow(target_chr_len)) {
		pre_len = sum(as.numeric(target_chr_len$len[0:(i-1)]))
		
		chr_lab_pos=pre_len + target_chr_len$len[i] %/% 2
		
		axis(3, at=chr_lab_pos, labels=target_chr_len$chr[i],tick=F, cex.axis=.5, line=-1)
		abline(v=pre_len, lty=3, col='gray')
		# lines( ( a$x[ which( a$space == chr[i] ) ] + pre_len ), a$copy[which(a$space==chr[i])], type="p", col = cols[as.numeric(as.character(a$state))] )
		
		points( ( a$x[ which( a$space == chr[i] ) ] + pre_len ), ( a$copy[which(a$space==chr[i])] + stat_offset ), col = cols[as.numeric(as.character(a$state[ which( a$space == chr[i] ) ]))], cex=.25)
		
		
		
		# lines(c(b$start[1], b$end[1]), rep(b$median[1], 2), lwd = 5, 
				# col = "green")
		# for (k in 1:nrow(b)) {
		
			# lines(c(b$start[k], b$end[k]), rep(b$median[k], 2), lwd = 5, 
				# col = "green")
		
		# }	
		
		# pre_len = sum(as.numeric(target_chr_len$len[0:(i-1)]))
		selected_chr=which( b$chr==target_chr_len$chr[i])
		medi=matrix(c(c(b$start[selected_chr], b$end[selected_chr])+pre_len, b$median[selected_chr]), ncol=3)
		
		# matrix(c( c(b$start[ which( b$chr==target_chr_len$chr[1] ) ], b$end[which(b$chr==target_chr_len$chr[1])] ) + pre_len), b$median[ which( b$chr==target_chr_len$chr[1] ) ])
		
		for (k in 1:nrow(medi)) {
			lines( c(medi[k,1], medi[k,2]), rep(medi[k,3]+ stat_offset, 2), col = "black", lwd=2 )
		}
	
	}
	
	
	# -------------------------------------
	chrs=chr
	# chrs='chr1'
	
	hs_cytoband=paste(Sys.getenv("data_path"), 'hs_cytoBand1.txt', sep='/')

	
cytoband0=read.table(hs_cytoband, col.names=c('chr', 'start', 'end', 'name', 'gie_stain'), stringsAsFactors=F)

	# cytoband=cytoband0[which(cytoband0$chr %in% chrs),] # if multi keys, use %in% instead of ==
	

	gie_stain=c("gneg",    "gpos25",  "gpos50",  "gpos75",  "gpos100", "acen",    "gvar",   "stalk")
	cols=c('white', 'gray75', 'gray50', 'gray25', 'black', 'red', 'cyan', 'blue');

	len=247249719

	y=c(0,1,1,0)

	y1=c(0,0.5,1)

	y2=c(0,0,0.5,1,1,0,0,0.5,1,1)


	# pdf("1.chrs.ideogram.pdf", width=12, height=2)

	# chr_num=cytoband$end[len0]/median(abs(cytoband$start-cytoband$end))

	# par('cin')

	# cin0=width*.8/chr_num*par('cin')/par('cin')[1]
	# par(cin=cin0)
	
	par(fig=c(0, 1, 0, 0.4), new=T)

	plot(xlim, c(0,1), type='n', ann=F, axes=F)
	
	# plot(NA, type='n', xlim=xlim, ylim=ylim1, cex.axis=.5)
	

	for (i in 1:nrow(target_chr_len)) {
		pre_len = sum(as.numeric(target_chr_len$len[0:(i-1)]))
		
		# chr_lab_pos=pre_len + target_chr_len$len[i] %/% 2
		cytoband=cytoband0[which(cytoband0$chr==target_chr_len$chr[i]),]
		len0=length(cytoband$chr)
		cytoband$start=cytoband$start+pre_len
		cytoband$end=cytoband$end+pre_len

		polygon(c(cytoband$start[1], cytoband$start[which(cytoband$gie_stain=='acen')[1]], cytoband$end[which(cytoband$gie_stain=='acen')[1]], cytoband$end[which(cytoband$gie_stain=='acen')[2]], cytoband$end[len0], cytoband$end[len0], cytoband$end[which(cytoband$gie_stain=='acen')[2]], cytoband$end[which(cytoband$gie_stain=='acen')[1]], cytoband$start[which(cytoband$gie_stain=='acen')[1]], cytoband$start[1]), y2)

		for ( i in 1:length(cytoband$chr)) {
			
			pq_pos=(cytoband$start[i] + cytoband$end[i]) %/% 2
			axis(1, at=pq_pos, labels=cytoband$name[i], tick=F, cex.axis=.2, line=-1, las=2)
			
			if ( cytoband$gie_stain[i]== "acen") {
				# polygon(c(x1,x1,x2,x2), y1, col=xx)
				# polygon(c(cytoband$start[i],cytoband$start[i],cytoband$end[i],cytoband$end[i]), y1, col=cols[which(gie_stain==cytoband$gie_stain[i])])
				polygon(c(cytoband$start[i],cytoband$end[i],cytoband$start[i]), y1, col=cols[which(gie_stain==cytoband$gie_stain[i])], border=NA)
				
				tmp=cytoband$start[i+1]
				cytoband$start[i+1]=cytoband$end[i+1]
				cytoband$end[i+1]=tmp
				
			}else {
				polygon(c(cytoband$start[i],cytoband$start[i],cytoband$end[i],cytoband$end[i]), y, col=cols[which(gie_stain==cytoband$gie_stain[i])], border=NA)
			}

			

		}
}
	
	
	
	dev.off()
}
